#----------------------------------------------------------#
#  Democratic Resilience (Aurel Croissant and Lars Lott)   #
#----------------------------------------------------------#

# Title: "Democratic Resilience in the Twenty-First Century. Search for an Analytical Framework and Explorative Analysis" #
# Authors: "Croissant, Aurel and Lott, Lars", Heidelberg University and FAU Erlangen-Nürnberg
# date: 2025-05-07
# journal: Political Studies
# DOI: 10.1177/00323217251345779
# written under "R version 4.5.0 (2025-04-11 ucrt)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# load packages
library(here)
library(tidyverse)
library(magrittr)
library(reshape2)
library(ggpubr)
library(gapminder)
library(countrycode)
library(ggrepel)
library(ggthemes)
library(lubridate)
library(estimatr)
library(texreg)
library(readxl)
library(haven)
library(scales)
library(zoo)
library(MCMCpack)
library(sandwich)
library(brglm)
library(GJRM)
library(lmtest)


#---------------------------------------------------#
#--------------------load data ---------------------#
#---------------------------------------------------#

# set you working directory here #

## load V-Dem data V14 ##

vdem_full_v14 <- readRDS(file.path("data/vdem14", "V-Dem-CY-Full+Others-v14.rds"))

## load Resilience Capacity data ##

ds.basic.subset.2000 <- readRDS("data/ds.basic.subset.2000_subset.rds")

## load ERT data v14 ##

ert_data <- readRDS("data/ert/ert_data_v14_uturns.rds")

summary(ert_data)

table(ert_data$uturn_ep_id)


bmr_data <- read.csv("data/bmr/bmr_data_extended.csv")

bmr_data <- bmr_data %>% dplyr::select(-X)

#### Data Merging ####

vdem_subset_v14 <- vdem_full_v14 %>%
  dplyr::select(country_id, year, COWcode, country_name, starts_with("v2x_polyarchy"), starts_with("v2x_libdem"),
                starts_with("v2x_regime"), v2caviol, v2caassemb, v2cagenmob, v2caconmob, v2cademmob, v2caautmob, v2cauni, 
                v2canuni, v2caprotac, v2cafres, v2cafexch, v2cainsaut, v2casurv, v2clacfree, e_regionpol_6C,
                v2xca_academ, e_pt_coup, e_regionpol, v2x_regime, v2xlg_legcon, v2x_jucon, e_gdppc, e_pop)

ert_data <- ert_data %>%
  dplyr::select(-c(country_text_id, country_name, v2x_regime:v2x_polyarchy_codehigh, change_edi, first_v2x_polyarchy))

vdem_subset_v14 <- vdem_subset_v14 %>%
  left_join(ert_data, by =c("country_id", "year"))

vdem_subset_v14 <- vdem_subset_v14 %>%
  left_join(ds.basic.subset.2000, by =c("country_id", "year"))

vdem_subset_v14 <- vdem_subset_v14 %>%
  left_join(bmr_data, by = c("COWcode" = "ccode", "year"))

vdem_subset_v14 <- vdem_subset_v14 %>%
  distinct(country_id, year, .keep_all = T)

#### Descriptive Tables ####
#### Resilience Capacity and Autocratization ####

vdem_subset_v14_tables <- vdem_subset_v14 %>%
  group_by(country_text_id) %>%
  mutate(v2x_polyarchy_lag = lag(v2x_polyarchy)) %>%
  mutate(aut_ep_outcome_agg = case_when(aut_ep_outcome_agg == 1 ~ "Democratic breakdown", 
                                        aut_ep_outcome_agg == 2 ~ "No democratic breakdown", 
                                        aut_ep_outcome_agg == 3 ~ "Regressed autocracy", 
                                        aut_ep_outcome_agg == 4 ~ "Outcome censored")) %>%
  mutate(aut_ep_turnaround_id = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_id,
                                       ifelse(!is.na(uturn_ep_id), uturn_ep_id, NA)), 
         aut_ep_turnaround_start_year = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_start_year,
                                       ifelse(!is.na(uturn_ep_id), uturn_ep_start, NA)), 
         aut_ep_turnaround_end_year = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_end_year,
                                               ifelse(!is.na(uturn_ep_id), uturn_ep_end, NA)), 
         aut_ep_turnaround_outcome = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_outcome_agg,
                                               ifelse(!is.na(uturn_ep_id), uturn_ep_outcome, NA))) 
table(vdem_subset_v14_tables$uturn_ep_id)
length(table(vdem_subset_v14_tables$uturn_ep_id)) # 108 episodes 

vdem_subset_v14_tables_df <- vdem_subset_v14_tables %>%
  filter(year >=2000) %>%
  filter(democracy == 1) %>%
  dplyr::select(country_text_id, country_name, year, starts_with("v2x_polyarchy"), v2x_regime, democracy, aut_ep_turnaround_id, aut_ep_turnaround_start_year, aut_ep_turnaround_end_year, 
                aut_ep_turnaround_outcome, ends_with("ResCap"), e_regionpol) %>%
  mutate(e_regionpol = case_when(e_regionpol == 1 ~"Eastern Europe and Central Asia", 
                                 e_regionpol == 2 ~"Latin America and the Caribbean", 
                                 e_regionpol == 3 ~"MENA", 
                                 e_regionpol == 4 ~"SSA", 
                                 e_regionpol == 5 ~"Western Europe and North America", 
                                 e_regionpol == 6 ~"Asia and Pacific"), 
         v2x_regime = case_when(v2x_regime == 1 ~ "Electoral Autocracy", 
                                  v2x_regime == 2 ~ "Electoral Democracy", 
                                  v2x_regime == 3 ~ "Liberal Democracy"))

writexl::write_xlsx(vdem_subset_v14_tables_df, "vdem_subset_v14_tables_df.xlsx")
length(table(vdem_subset_v14_tables_df$aut_ep_turnaround_id)) # 68 

country_names <- vdem_subset_v14_tables_df %>%
         drop_na(ResCap) %>%
          ungroup() %>%
         distinct(country_name)# 103 

table(vdem_subset_v14_tables_df$aut_ep_turnaround_id)

## List of all Autocratization Episodes ##

table_5 <- vdem_subset_v14_tables %>%
  drop_na(ResCap) %>%
  filter(!is.na(aut_ep_turnaround_id)) %>%
  group_by(aut_ep_turnaround_id) %>%
  summarize(country = first(country_name), 
            start_year = first(aut_ep_turnaround_start_year), 
            end_year = first(aut_ep_turnaround_end_year), 
            outcome = last(aut_ep_turnaround_outcome)) 

table_5$start_year <- as.character(table_5$start_year)
table_5$end_year <- as.character(table_5$end_year)

table_5 <- table_5 %>%
  arrange(outcome, start_year)

modelsummary::datasummary_df(data = table_5, output = '../results/Table_A2.docx', fmt= 3)

table(table_5$outcome)

table_5 <- table_5 %>%
  arrange(country, start_year)

#### Table 7: Resilience Capacity 

table5 <- vdem_subset_v14_tables %>%
  group_by(country_text_id) %>%
  mutate(ResCap_lag = lag(ResCap)) %>%
  drop_na(ResCap) %>%
  group_by(aut_ep_turnaround_id) %>%
  fill(aut_ep_turnaround_outcome, .direction = "up") %>%
  summarize(ResCap_lag = first(ResCap), 
            aut_ep_turnaround_outcome = first(aut_ep_turnaround_outcome)) %>%
  group_by(aut_ep_turnaround_outcome) %>%
  summarize(mean = mean(ResCap_lag, na.rm = T), 
            sd = sd(ResCap_lag, na.rm = T), 
            min = min(ResCap_lag, na.rm = T), 
            max = max(ResCap_lag, na.rm = T)) %>%
  drop_na(aut_ep_turnaround_outcome)
modelsummary::datasummary_df(data = table5, output = '../results/Table_4.docx', fmt= 3)

table3 <- vdem_subset_v14_tables %>%
  group_by(country_text_id) %>%
  mutate(ResCap_lag = lag(ResCap)) %>%
  drop_na(ResCap) %>%
  group_by(aut_ep_turnaround_id) %>%
  fill(aut_ep_turnaround_outcome, .direction = "up") %>%
  summarize(ResCap_lag = first(ResCap), 
            aut_ep_turnaround_outcome = first(aut_ep_turnaround_outcome)) %>%
  group_by(aut_ep_turnaround_outcome) %>%
  summarize( n = as.integer(n()))

modelsummary::datasummary_df(data = table3, output = '../results/Table_3.docx', fmt= 0)



#################################################################################
#################################################################################
#################################################################################

#### Data Wrangling ####

vdem_subset_v14 <- vdem_subset_v14 %>%
  group_by(e_regionpol, year) %>%
  mutate(regionpol_EDI = mean(v2x_polyarchy, na.rm = TRUE)) %>%
  ungroup() %>%
  group_by(country_id) %>%
  mutate(e_gdppc_growth = ((e_gdppc - dplyr::lag(e_gdppc))/dplyr::lag(e_gdppc))*100) %>%
  mutate(e_gdppc_growth_roll5 = rollmean(e_gdppc_growth, k = 5, fill = NA)) %>%
  ungroup() %>%
  group_by(aut_ep_id) %>%
  mutate(aut_duration =  year- first(year), 
         aut_duration = ifelse(is.na(aut_ep_id), 0, aut_duration)) %>%
  group_by(country_id) %>%
  mutate(group_id = cumsum(aut_ep != lag(aut_ep, default = FALSE))) %>% # group_id for each autocra_period
  group_by(country_id, group_id) %>%
  mutate(spell = year-first(year), 
         spell = ifelse(is.na(aut_ep_id), spell, 0))

vdem_subset_v14 <- vdem_subset_v14 %>%
  distinct(country_text_id, year, .keep_all = T)


#---------------------------------------------------#
#------------Democracy Stock Variable---------------#
#---------------------------------------------------#

decay_sum <- function(tm, vl, decay) {
  last_time <- 0
  current_sum <- 0
  sums <- numeric(length(vl))
  ldecay <- (1-decay)
  for (i in 1:length(vl)) {
    delta <- as.numeric(tm[i] - last_time)
    current_sum <- current_sum * (ldecay * delta) + vl[i]
    last_time <- tm[i]
    sums[i] <- current_sum
  }
  sums
}

annual_depricition_function  <- function(year, var, decay) {
  last_time <- first(year)
  current_sum <- 0
  sums <- numeric(length(var))
  ldecay <- (1-decay)
  for (i in 1:length(var)) {
    delta <- as.numeric(year[i] - last_time)
    current_sum <- ldecay * (current_sum + decay * var[i])
    last_time <- year[i]
    sums[i] <- current_sum
  }
  sums
}


## Functions ##

na_interpolation2 <- function(x, option = "linear", maxgap) {
  library(imputeTS)
  library(dplyr)
  
  total_not_missing <- sum(!is.na(x))
  
  # check there is sufficient data for na_interpolation 
  if(total_not_missing < 2) {x} 
  
  else
    
    # replace takes an input vector, a T/F vector & replacement value
  {replace(
    # input vector is interpolated data
    # this will impute leading/lagging NAs which we don't want 
    imputeTS::na_interpolation(x, option = option, maxgap = maxgap), 
    
    # create T/F vector for NAs,  
    is.na(na.approx(x, na.rm = FALSE)), 
    
    # replace TRUE with NA in input vector  
    NA) 
  }
}

vdem_subset_v14 <- vdem_subset_v14 %>%
  group_by(country_id) %>%
  mutate(v2x_polyarchy_interpolated = na_interpolation2(v2x_polyarchy, option = "linear", maxgap = 5)) %>%
  drop_na(v2x_polyarchy_interpolated) %>%
  mutate(democracy_stock = annual_depricition_function(year, v2x_polyarchy_interpolated, 0.01)) %>%
  dplyr::select(country_id, country_name, year, v2x_polyarchy, democracy_stock, everything())

summary(vdem_subset_v14$e_pop)
summary(vdem_subset_v14$e_gdppc)

vdem_subset_v14 <- vdem_subset_v14 %>%
  group_by(country_id) %>%
  fill(e_pop, e_gdppc, e_gdppc_growth, e_gdppc_growth_roll5)


#### Breakdown Resilience (Second Stage) #####

vdem_subset_v14 <- vdem_subset_v14 %>%
  group_by(country_text_id) %>%
  mutate(v2x_regime_dem = ifelse(v2x_regime >=2, 1, 0), 
         dem_breakdown_v2x_regime = ifelse((lag(v2x_regime_dem) - v2x_regime_dem) == 1, 1, 0)) %>%
  mutate(stage2_res = ifelse(aut_ep == 1 & dem_breakdown_v2x_regime == 1, 1, 0))

test <- vdem_subset_v14 %>%
  dplyr::select(country_text_id, year, aut_ep, stage2_res, v2x_regime_dem, dem_breakdown_v2x_regime) %>%
  filter(year >=1900)



## lags for independent variables ##

vdem_subset_v14_noNA <- vdem_subset_v14 %>%
  mutate(e_gdppclog = log(e_gdppc+1), 
         e_poplog = log(e_pop +1), 
         year_1900 = year -1900) %>%
  group_by(country_id) %>%
  mutate(v2cademmob_lag = lag(v2cademmob), 
         e_gdppclog_lag = lag(e_gdppclog), 
         e_gdppc_growth_roll5_lag = lag(e_gdppc_growth_roll5), 
         regionpol_EDI_lag = lag(regionpol_EDI), 
         e_poplog_lag = lag(e_poplog), 
         v2x_polyarchy_lag = lag(v2x_polyarchy), 
         v2x_polyarchy_lag_squared = v2x_polyarchy_lag^2, 
         v2x_jucon_lag = lag(v2x_jucon), 
         v2xlg_legcon_lag = lag(v2xlg_legcon), 
         democracy_stock_lag = lag(democracy_stock), 
         ResCap_lag = lag(ResCap), 
         AResCap_lag = lag(AResCap), 
         MResCap_lag = lag(MResCap))

vdem_subset_v14_noNA <- vdem_subset_v14_noNA %>%
  mutate(region_factor = as.factor(e_regionpol_6C),
         region_factor = relevel(region_factor, ref=1), 
         spell = spell/100,
         spell_sq = spell^2/100,
         spell_cb = spell^3/100,
         year_sq = year_1900^2)


## drop NA ## 
vdem_subset_v14_noNA <- vdem_subset_v14_noNA %>%
  group_by(country_text_id) %>%
  mutate(v2x_regime_lag = lag(v2x_regime)) %>%
  ungroup() %>%
  filter(year >=2000) %>%
  filter(democracy == 1) %>%
  mutate(aut_ep_stage2 = ifelse(aut_ep_outcome %in% c(1,5), 1, 0)) %>%
  drop_na(ResCap_lag, e_gdppclog, e_gdppc_growth_roll5_lag, regionpol_EDI_lag, 
          e_poplog,region_factor ,v2x_jucon_lag, v2xlg_legcon_lag, year_1900, spell, country_id, democracy_stock_lag) 

summary(vdem_subset_v14_noNA)
length(unique(vdem_subset_v14_noNA$country_name)) # 103 countries 

## List of all Autocratization Episodes ##

vdem_subset_v14_noNA <- vdem_subset_v14_noNA %>%
  group_by(country_text_id) %>%
  mutate(v2x_polyarchy_lag = lag(v2x_polyarchy)) %>%
  mutate(aut_ep_outcome_agg = case_when(aut_ep_outcome_agg == 1 ~ "Democratic breakdown", 
                                        aut_ep_outcome_agg == 2 ~ "No democratic breakdown", 
                                        aut_ep_outcome_agg == 3 ~ "Regressed autocracy", 
                                        aut_ep_outcome_agg == 4 ~ "Outcome censored")) %>%
  mutate(aut_ep_turnaround_id = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_id,
                                       ifelse(!is.na(uturn_ep_id), uturn_ep_id, NA)), 
         aut_ep_turnaround_start_year = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_start_year,
                                               ifelse(!is.na(uturn_ep_id), uturn_ep_start, NA)), 
         aut_ep_turnaround_end_year = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_end_year,
                                             ifelse(!is.na(uturn_ep_id), uturn_ep_end, NA)), 
         aut_ep_turnaround_outcome = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_outcome_agg,
                                            ifelse(!is.na(uturn_ep_id), uturn_ep_outcome, NA)))

table_resCap_desc <- vdem_subset_v14_noNA %>%
  group_by(aut_ep_turnaround_id) %>%
  summarize(country = first(country_name), 
            start_year = first(aut_ep_turnaround_start_year), 
            end_year = first(aut_ep_turnaround_end_year), 
            outcome = last(aut_ep_turnaround_outcome), 
            ResCap = first(ResCap_lag), 
            region =first(region_factor), 
            start_regime = first(v2x_regime_lag)) %>%
  drop_na(aut_ep_turnaround_id)

table_resCap_desc$start_year <- as.character(table_resCap_desc$start_year)
table_resCap_desc$end_year <- as.character(table_resCap_desc$end_year)

table_resCap_desc <- table_resCap_desc %>%
  mutate(region = case_when(region == 1 ~"Eastern Europe and Central Asia", 
                            region == 2 ~"Latin America and the Caribbean", 
                            region == 3 ~"MENA", 
                            region == 4 ~"SSA", 
                            region == 5 ~"Western Europe and North America", 
                            region == 6 ~"Asia and Pacific"), 
         start_regime = case_when(start_regime == 1 ~ "Electoral Autocracy", 
                                  start_regime == 2 ~ "Electoral Democracy", 
                                  start_regime == 3 ~ "Liberal Democracy"))

table_resCap_desc <- table_resCap_desc %>%
  arrange(outcome, start_year)

modelsummary::datasummary_df(data = table_resCap_desc, output = '../results/Table_Des_ResCap.docx', fmt= 3)

## setting ongoing events to NA as recommend by McGrath (2015)

vdem_subset_v14_noNA_stage2 <- vdem_subset_v14_noNA %>%
  group_by(country_id) %>%
  mutate(aut_ep_onset = ifelse(aut_ep==1 & lag(aut_ep)== 1, NA, aut_ep)) %>%
  dplyr::select(country_id, year, aut_ep, aut_ep_onset, everything())

vdem_subset_v14_noNA <- vdem_subset_v14_noNA %>%
  group_by(country_id) %>%
  mutate(aut_ep_onset = ifelse(aut_ep==1 & lag(aut_ep)== 1, NA, aut_ep)) %>%
  dplyr::select(country_id, year, aut_ep, aut_ep_onset, everything()) %>%
  drop_na(aut_ep_onset)

summary(vdem_subset_v14_noNA$aut_ep_onset) # 3.03 % of all country-years
table(vdem_subset_v14_noNA$aut_ep_onset) # 47 autocratization onsets after 2000

vdem_subset_v14_noNA <- vdem_subset_v14_noNA %>%
  mutate(stage1_res = ifelse(aut_ep == 0, 1, 0))

summary(vdem_subset_v14_noNA$stage1_res) # 96.98% of all country-years

# list of all onsets #

list_of_auto_onsets <- vdem_subset_v14_noNA %>%
  filter(aut_ep_onset == 1) %>%
  ungroup() %>%
  dplyr::select(country_name, year)

list_of_auto_onsets$year <- as.factor(list_of_auto_onsets$year)

modelsummary::datasummary_df(data = list_of_auto_onsets, 
               output = "../results/TableA2.docx", 
               title = "List of country-years with autocratization onset", 
               longtable=TRUE)

country_list_auto <- unique(list_of_auto_onsets$country_name)


# list of all onsets with breakdown  #

list_of_auto_onsets <- vdem_subset_v14_noNA %>%
  filter(aut_ep_onset == 1) %>%
  ungroup() %>%
  dplyr::select(country_name, year, aut_ep_start_year, aut_ep_end_year, aut_ep_outcome) %>%
  mutate(aut_ep_outcome = case_when(aut_ep_outcome == 1 ~ "Democratic Breakdown", 
                                      aut_ep_outcome == 2 ~ "Preempted Democratic Breakdown", 
                                      aut_ep_outcome == 3 ~ "Diminished Democracy", 
                                      aut_ep_outcome == 4 ~ "Averted Regression",
                                      aut_ep_outcome == 5 ~ "Regressed Autocracy",
                                      aut_ep_outcome == 6 ~ "Outcome Censored"))

list_of_auto_onsets$year <- as.factor(list_of_auto_onsets$year)
list_of_auto_onsets$aut_ep_start_year <- as.factor(list_of_auto_onsets$aut_ep_start_year)
list_of_auto_onsets$aut_ep_end_year <- as.factor(list_of_auto_onsets$aut_ep_end_year)

modelsummary::datasummary_df(data = list_of_auto_onsets, 
                             output = "../results/TableA3.docx", 
                             title = "List of country-years with autocratization onset", 
                             longtable=TRUE)

country_list_auto <- unique(list_of_auto_onsets$country_name)

vdem_subset_v14_noNA_resap <- vdem_subset_v14_noNA %>%
  dplyr::select(country_text_id, country_name, year, aut_ep_onset, ResCap)

##################################################################################
##################################################################################

#### Linear Models #####

# Run model 1
onset_m1 <- brglm::brglm(as.factor(aut_ep_onset) ~ ResCap_lag + 
                           e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                           regionpol_EDI_lag + 
                           region_factor + year_1900 + year_sq, 
                         data = vdem_subset_v14_noNA, family = binomial(link = "probit"),
                         method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                         pl=T)
summary(onset_m1)

#Clustered SEs
ses_ons1_cl <-  coeftest(onset_m1, vcov = vcovCL(onset_m1, cluster = vdem_subset_v14_noNA[, "country_id"]),
                         type = "HC1", fix = T, cadjust = T)
#ses_ons1_cl <- tidy(ses_ons1_cl) 


# Run model 2
onset_m2 <- brglm::brglm(as.factor(aut_ep_onset) ~ AResCap_lag + 
                           e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                           regionpol_EDI_lag + 
                           region_factor + year_1900 + year_sq, 
                         data = vdem_subset_v14_noNA, family = binomial(link = "probit"),
                         method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                         pl=T)
summary(onset_m2)

#Clustered SEs
ses_ons2_cl <-  coeftest(onset_m2, vcov = vcovCL(onset_m2, vdem_subset_v14_noNA[, "country_id"]),
                         type = "HC1", fix = T, cadjust = T)
#ses_ons2_cl <- tidy(ses_ons2_cl) 


# Run model 3
onset_m3 <- brglm::brglm(as.factor(aut_ep_onset) ~ MResCap_lag +
                           e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                           regionpol_EDI_lag + 
                           region_factor + year_1900 + year_sq, 
                         data = vdem_subset_v14_noNA, family = binomial(link = "probit"),
                         method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                         pl=T)
summary(onset_m3)

#Clustered SEs
ses_ons3_cl <-  coeftest(onset_m3, vcov = vcovCL(onset_m3, vdem_subset_v14_noNA[, "country_id"]),
                         type = "HC1", fix = T, cadjust = T)
#ses_ons3_cl <- tidy(ses_ons3_cl) 

## Extract models ##
texreg::wordreg(list(onset_m1, onset_m2, onset_m3), 
                override.pvalues=list(ses_ons1_cl[, 4], ses_ons2_cl[, 4], ses_ons3_cl[, 4]),
                override.se= list(ses_ons1_cl[, 2], ses_ons2_cl[, 2], ses_ons3_cl[, 2]), 
                file = "../results/Table_A3.docx", 
                custom.coef.names=c("Intercept", 
                                    "Democratic Resilience Capacity", 
                                    "GDP pc log",
                                    "GDP growth",
                                    "Population log",
                                    "Regional democracy levels",
                                    "Western Europe and North America",
                                    "Subsaharan Africa",
                                    "Asia and Pacific", 
                                    "Eastern Europe and Central Asia",
                                    "MENA", 
                                    "Year",
                                    "Year squared",
                                    "Additive DRC", 
                                    "Mulitplicative DRC"))

# Plot marginal effects 

af_onset <-  sjPlot::plot_model(onset_m1, type = "eff",
                                terms = c("ResCap_lag[all]"),
                                vcov.fun = "vcovCL", 
                                vcov.type = "HC1",
                                vcov.args = "cluster = vdem_subset_v14_noNA$country_id",
                                alpha = 0.3) + 
  ylab("P(autocratization onset)") + xlab("DRC Index") + theme_bw() + 
  scale_y_continuous(limits = c(0,0.2), breaks = seq(0,0.2,0.05)) +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0,1,0.2)) + ggtitle("") +
  geom_rug(data = vdem_subset_v14_noNA, aes(x = ResCap_lag, y = 0), alpha = .1, position = "jitter", sides = "b")

af_onset_data <- af_onset$data

data_description <- vdem_subset_v14_noNA %>% 
 # dplyr::filter(ResCap %in% c(0.19:0.3)) %>%
  dplyr::select(country_name, year, aut_ep, ResCap)
  


af_onset2 <-  sjPlot::plot_model(onset_m2, type = "eff",
                                 terms = c("AResCap_lag[all]"),
                                 vcov.fun = "vcovCL", 
                                 vcov.type = "HC1",
                                 vcov.args = "cluster = vdem_subset_v14_noNA$country_id",
                                 alpha = 0.3) + 
  ylab("P(autocratization onset)") + xlab("Additive DRC Index") + theme_bw() + 
  scale_y_continuous(limits = c(0,0.2), breaks = seq(0,0.2,0.05)) +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0,1,0.2)) + ggtitle("") +
  geom_rug(data = vdem_subset_v14_noNA, aes(x = AResCap_lag, y = 0), alpha = .1, position = "jitter", sides = "b")


af_onset3 <-  sjPlot::plot_model(onset_m3, type = "eff",
                                 terms = c("MResCap_lag[all]"),
                                 vcov.fun = "vcovCL", 
                                 vcov.type = "HC1",
                                 vcov.args = "cluster = vdem_subset_v14_noNA$country_id",
                                 alpha = 0.3) + 
  ylab("P(autocratization onset)") + xlab("Multiplicative DRC Index") + theme_bw() + 
  scale_y_continuous(limits = c(0,0.2), breaks = seq(0,0.2,0.05)) +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0,1,0.2)) + ggtitle("") +
  geom_rug(data = vdem_subset_v14_noNA, aes(x = MResCap_lag, y = 0), alpha = .1, position = "jitter", sides = "b")

figure <- ggpubr::ggarrange(af_onset, af_onset2,af_onset3,
                            nrow = 2, ncol = 2, legend = "bottom", 
                            labels = c("A","B", "C"))


ggsave("../results/Figure_5.png", width = 30,
       height = 30,
       units = c("cm"), dpi = 1000)

ggsave("../results/Figure_5_low_res.png", width = 30,
       height = 30,
       units = c("cm"), dpi = 400)


##################################################################################
##################################################################################

#### Dimensions of Resilience Capacity #####

#### Data Analysis ####
## onset model ##
# Rare events logistic regression using the brglm package

# Run model 5: Inst Var 
onset_m4 <- brglm::brglm(as.factor(aut_ep_onset) ~ inst_dim + #democracy_stock_lag +
                           e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                           regionpol_EDI_lag + 
                           #v2x_jucon_lag + v2xlg_legcon_lag +
                           region_factor + year_1900 + year_sq, 
                         data = vdem_subset_v14_noNA, family = binomial(link = "probit"),
                         method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                         pl=T)
summary(onset_m4)

#Clustered SEs
ses_ons4_cl <-  coeftest(onset_m4, vcov = vcovCL(onset_m4, cluster = vdem_subset_v14_noNA[, "country_id"]),
                         type = "HC1", fix = T, cadjust = T)


# Run model 5: Pol Party Var
onset_m5 <- brglm::brglm(as.factor(aut_ep_onset) ~ polpar_dim + #democracy_stock_lag +
                           e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                           regionpol_EDI_lag + 
                           #v2x_jucon_lag + v2xlg_legcon_lag +
                           region_factor + year_1900 + year_sq, 
                         data = vdem_subset_v14_noNA, family = binomial(link = "probit"),
                         method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                         pl=T)
summary(onset_m5)

#Clustered SEs
ses_ons5_cl <-  coeftest(onset_m5, vcov = vcovCL(onset_m5, vdem_subset_v14_noNA[, "country_id"]),
                         type = "HC1", fix = T, cadjust = T)


# Run model 6: Societal Var
onset_m6 <- brglm::brglm(as.factor(aut_ep_onset) ~ civsoc_dim + #democracy_stock_lag +
                           e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                           regionpol_EDI_lag + 
                           #v2x_jucon_lag + v2xlg_legcon_lag +
                           region_factor + year_1900 + year_sq, 
                         data = vdem_subset_v14_noNA, family = binomial(link = "probit"),
                         method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                         pl=T)
summary(onset_m6)

#Clustered SEs
ses_ons6_cl <-  coeftest(onset_m6, vcov = vcovCL(onset_m6, vdem_subset_v14_noNA[, "country_id"]),
                         type = "HC1", fix = T, cadjust = T)

# Run model 7: Societal Var
onset_m7 <- brglm::brglm(as.factor(aut_ep_onset) ~ polcom_dim + #democracy_stock_lag +
                           e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                           regionpol_EDI_lag + 
                           #v2x_jucon_lag + v2xlg_legcon_lag +
                           region_factor + year_1900 + year_sq, 
                         data = vdem_subset_v14_noNA, family = binomial(link = "probit"),
                         method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                         pl=T)
summary(onset_m7)

#Clustered SEs
ses_ons7_cl <-  coeftest(onset_m7, vcov = vcovCL(onset_m7, vdem_subset_v14_noNA[, "country_id"]),
                         type = "HC1", fix = T, cadjust = T)

# Run model 8: All Dimensions 
onset_m8 <- brglm::brglm(as.factor(aut_ep_onset) ~ inst_dim + 
                           polpar_dim + civsoc_dim + polcom_dim +
                           e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                           regionpol_EDI_lag + 
                           #v2x_jucon_lag + v2xlg_legcon_lag +
                           region_factor + year_1900 + year_sq, 
                         data = vdem_subset_v14_noNA, family = binomial(link = "probit"),
                         method = "brglm.fit", control.brglm = brglm.control(br.maxit = 1000, br.trace = T),
                         pl=T)
summary(onset_m8)

#Clustered SEs
ses_ons8_cl <-  coeftest(onset_m8, vcov = vcovCL(onset_m8, vdem_subset_v14_noNA[, "country_id"]),
                         type = "HC1", fix = T, cadjust = T)


## Extract models ##
texreg::wordreg(list(onset_m4, onset_m5, onset_m6, onset_m7, onset_m8), 
                override.pvalues=list(ses_ons4_cl[, 4], ses_ons5_cl[, 4], ses_ons6_cl[, 4], ses_ons7_cl[, 4], ses_ons8_cl[, 4]),
                override.se= list(ses_ons4_cl[, 2], ses_ons5_cl[, 2], ses_ons6_cl[, 2], ses_ons7_cl[, 2], ses_ons8_cl[, 2]), 
                file = "../results/Table_A5.docx", 
                custom.coef.names=c("Intercept", 
                                    "Institutional Dimension", 
                                    "GDP pc log",
                                    "GDP growth",
                                    "Population log",
                                    "Regional democracy levels",
                                    "Western Europe and North America",
                                    "Subsaharan Africa",
                                    "Asia and Pacific", 
                                    "Eastern Europe and Central Asia",
                                    "MENA", 
                                    "Year",
                                    "Year squared",
                                    "Political Parties Dimension", 
                                    "Civil society and civic culture Dimension", 
                                    "Political Community Dimension"))

# Plot marginal effects 

m4_pred <-  sjPlot::plot_model(onset_m4, type = "eff",
                                terms = c("inst_dim[all]"),
                                vcov.fun = "vcovCL", 
                                vcov.type = "HC1",
                                vcov.args = "cluster = vdem_subset_v14_noNA$country_id",
                                alpha = 0.3) + 
  ylab("P(autocratization onset)") + xlab("Institutional Dimension") + theme_bw() + 
  scale_y_continuous(limits = c(0,0.2), breaks = seq(0,0.2,0.05)) +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0,1,0.2)) + ggtitle("") +
  geom_rug(data = vdem_subset_v14_noNA, aes(x = inst_dim, y = 0), alpha = .1, position = "jitter", sides = "b")


m5_pred <-  sjPlot::plot_model(onset_m5, type = "eff",
                                 terms = c("polpar_dim[all]"),
                                 vcov.fun = "vcovCL", 
                                 vcov.type = "HC1",
                                 vcov.args = "cluster = vdem_subset_v14_noNA$country_id",
                                 alpha = 0.3) + 
  ylab("P(autocratization onset)") + xlab("Political Parties Dimension") + theme_bw() + 
  scale_y_continuous(limits = c(0,0.2), breaks = seq(0,0.2,0.05)) +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0,1,0.2)) + ggtitle("") +
  geom_rug(data = vdem_subset_v14_noNA, aes(x = polpar_dim, y = 0), alpha = .1, position = "jitter", sides = "b")


m6_pred <-  sjPlot::plot_model(onset_m6, type = "eff",
                                 terms = c("civsoc_dim[all]"),
                                 vcov.fun = "vcovCL", 
                                 vcov.type = "HC1",
                                 vcov.args = "cluster = vdem_subset_v14_noNA$country_id",
                                 alpha = 0.3) + 
  ylab("P(autocratization onset)") + xlab("Civil society and civic culture Dimension") + theme_bw() + 
  scale_y_continuous(limits = c(0,0.2), breaks = seq(0,0.2,0.05)) +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0,1,0.2)) + ggtitle("") +
  geom_rug(data = vdem_subset_v14_noNA, aes(x = civsoc_dim, y = 0), alpha = .1, position = "jitter", sides = "b")

m7_pred <-  sjPlot::plot_model(onset_m7, type = "eff",
                               terms = c("polcom_dim[all]"),
                               vcov.fun = "vcovCL", 
                               vcov.type = "HC1",
                               vcov.args = "cluster = vdem_subset_v14_noNA$country_id",
                               alpha = 0.3) + 
  ylab("P(autocratization onset)") + xlab("CPolitical Community Dimension") + theme_bw() + 
  scale_y_continuous(limits = c(0,0.2), breaks = seq(0,0.2,0.05)) +
  scale_x_continuous(limits = c(0, 1), breaks = seq(0,1,0.2)) + ggtitle("") +
  geom_rug(data = vdem_subset_v14_noNA, aes(x = polcom_dim, y = 0), alpha = .1, position = "jitter", sides = "b")

figure <- ggpubr::ggarrange(m4_pred, m5_pred, m6_pred, m7_pred, 
                            nrow = 2, ncol = 2, legend = "bottom", 
                            labels = c("A","B", "C", "D"))


ggsave("../results/Figure_A12.png", width = 30,
       height = 30,
       units = c("cm"), dpi = 1000)

ggsave("../results/Figure_A12_low_res.png", width = 30,
       height = 30,
       units = c("cm"), dpi = 400)

################################################################################
################################################################################

#### Table A4 ####
summary(vdem_subset_v14_noNA_stage2$stage2_res) # 0.9% of all country-years resulted in a democratic breakdown

# run selection model using gjrm package

selection_m1 <- gjrm(list(as.factor(aut_ep) ~ ResCap_lag  +
                            e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                            regionpol_EDI_lag +
                            region_factor + year_1900 + year_sq, 
                          as.factor(stage2_res) ~ ResCap_lag  +
                            e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                            regionpol_EDI_lag +
                            region_factor + year_1900 + year_sq), 
                     data = vdem_subset_v14_noNA_stage2, model = "BSS", copula = "C0",
                     margins = c("probit", "probit"))
summary(selection_m1)


data_description <- vdem_subset_v14_noNA_stage2 %>% 
  # dplyr::filter(ResCap %in% c(0.19:0.3)) %>%
  dplyr::select(country_name, year, aut_ep, stage2_res, ResCap)

# run selection model using gjrm package

selection_m2 <- gjrm(list(as.factor(aut_ep) ~ MResCap_lag  +
                            e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                            regionpol_EDI_lag +
                            region_factor + year_1900 + year_sq, 
                          as.factor(stage2_res) ~ MResCap_lag  +
                            e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                            regionpol_EDI_lag +
                            region_factor + year_1900 + year_sq), 
                     data = vdem_subset_v14_noNA_stage2, model = "BSS", copula = "C0",
                     margins = c("probit", "probit"))
summary(selection_m2)


# run selection model using gjrm package

selection_m3 <- gjrm(list(as.factor(aut_ep) ~ AResCap_lag  +
                            e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                            regionpol_EDI_lag +
                            region_factor + year_1900 + year_sq, 
                          as.factor(stage2_res) ~ AResCap_lag  +
                            e_gdppclog_lag + e_gdppc_growth_roll5_lag +  e_poplog_lag +
                            regionpol_EDI_lag +
                            region_factor + year_1900 + year_sq), 
                     data = vdem_subset_v14_noNA_stage2, model = "BSS", copula = "C0",
                     margins = c("probit", "probit"))
summary(selection_m3)

conv.check(selection_m1)    
conv.check(selection_m2)                          
conv.check(selection_m3)                          

mean(selection_m1$theta)
mean(selection_m2$theta)
mean(selection_m3$theta)

selection_m1 <- adjCov(selection_m1, id = vdem_subset_v14_noNA_stage2$country_text_id)
selection_m2 <- adjCov(selection_m2, id = vdem_subset_v14_noNA_stage2$country_text_id)
selection_m3 <- adjCov(selection_m3, id = vdem_subset_v14_noNA_stage2$country_text_id)

summary(selection_m1)

##### Extracting Table 8: Main results #####

## M1 Selection Model ##

sumamry_selection_m1 <- summary(selection_m1)

model2_1st_coef <- sumamry_selection_m1$tableP1[,1]
model2_2nd_coef <- sumamry_selection_m1$tableP2[,1]

model2_1st_se <- sumamry_selection_m1$tableP1[,2]
model2_2nd_se <- sumamry_selection_m1$tableP2[,2]

model2_1st_p <- sumamry_selection_m1$tableP1[,4]
model2_2nd_p <- sumamry_selection_m1$tableP2[,4]

# Model 2: First Stage 
model2_1st_coef <- bind_rows(lapply(model2_1st_coef, as.data.frame), .id = "model")
model2_1st_se <- bind_rows(lapply(model2_1st_se, as.data.frame), .id = "model")
model2_1st_p <- bind_rows(lapply(model2_1st_p, as.data.frame), .id = "model")

names(model2_1st_coef) <- c( "term", "coef")
names(model2_1st_se) <- c( "term", "se")
names(model2_1st_p) <- c( "term", "pval")

model2_1st <- model2_1st_coef %>%
  full_join(model2_1st_se, by = c("term")) %>%
  full_join(model2_1st_p, by = c("term"))

# Model 2: Second Stage
model2_2nd_coef <- bind_rows(lapply(model2_2nd_coef, as.data.frame), .id = "model")
model2_2nd_se <- bind_rows(lapply(model2_2nd_se, as.data.frame), .id = "model")
model2_2nd_p <- bind_rows(lapply(model2_2nd_p, as.data.frame), .id = "model")

names(model2_2nd_coef) <- c( "term", "coef")
names(model2_2nd_se) <- c( "term", "se")
names(model2_2nd_p) <- c( "term", "pval")

model2_2nd <- model2_2nd_coef %>%
  full_join(model2_2nd_se, by = c("term")) %>%
  full_join(model2_2nd_p, by = c("term"))

# Function to add stars based on p-values
add_stars <- function(pvalue) {
  ifelse(pvalue < 0.01, "***", 
         ifelse(pvalue < 0.05, "**", 
                ifelse(pvalue < 0.1, "*", "")))
}

model2_1st <- model2_1st %>%
  mutate(
    coef_stars = paste0(round(coef, 3), add_stars(pval)),
    se_paren = paste0("(", round(se, 3), ")")
  ) %>%
  pivot_longer(cols = c(coef_stars, se_paren), names_to = "type", values_to = "estimate") %>%
  arrange(term, type) %>%
  dplyr::select(term, estimate)

model2_2nd <- model2_2nd %>%
  mutate(
    coef_stars = paste0(round(coef, 3), add_stars(pval)),
    se_paren = paste0("(", round(se, 3), ")")
  ) %>%
  pivot_longer(cols = c(coef_stars, se_paren), names_to = "type", values_to = "estimate") %>%
  arrange(term, type) %>%
  dplyr::select(term, estimate)

model_second_stage <- cbind(model2_1st, model2_2nd[,2])

modelsummary::datasummary_df(data = model_second_stage, 
                             output = "../results/Table_A4_M2_1st.docx", 
                             title = "Model 2", longtable = T)


## M2 Selection Model ##

sumamry_selection_m2 <- summary(selection_m2)

model3_1st_coef <- sumamry_selection_m2$tableP1[,1]
model3_2nd_coef <- sumamry_selection_m2$tableP2[,1]

model3_1st_se <- sumamry_selection_m2$tableP1[,2]
model3_2nd_se <- sumamry_selection_m2$tableP2[,2]

model3_1st_p <- sumamry_selection_m2$tableP1[,4]
model3_2nd_p <- sumamry_selection_m2$tableP2[,4]

# Model 2: First Stage 
model3_1st_coef <- bind_rows(lapply(model3_1st_coef, as.data.frame), .id = "model")
model3_1st_se <- bind_rows(lapply(model3_1st_se, as.data.frame), .id = "model")
model3_1st_p <- bind_rows(lapply(model3_1st_p, as.data.frame), .id = "model")

names(model3_1st_coef) <- c( "term", "coef")
names(model3_1st_se) <- c( "term", "se")
names(model3_1st_p) <- c( "term", "pval")

model3_1st <- model3_1st_coef %>%
  full_join(model3_1st_se, by = c("term")) %>%
  full_join(model3_1st_p, by = c("term"))

# Model 2: Second Stage
model3_2nd_coef <- bind_rows(lapply(model3_2nd_coef, as.data.frame), .id = "model")
model3_2nd_se <- bind_rows(lapply(model3_2nd_se, as.data.frame), .id = "model")
model3_2nd_p <- bind_rows(lapply(model3_2nd_p, as.data.frame), .id = "model")

names(model3_2nd_coef) <- c( "term", "coef")
names(model3_2nd_se) <- c( "term", "se")
names(model3_2nd_p) <- c( "term", "pval")

model3_2nd <- model3_2nd_coef %>%
  full_join(model3_2nd_se, by = c("term")) %>%
  full_join(model3_2nd_p, by = c("term"))

# Function to add stars based on p-values
add_stars <- function(pvalue) {
  ifelse(pvalue < 0.01, "***", 
         ifelse(pvalue < 0.05, "**", 
                ifelse(pvalue < 0.1, "*", "")))
}

model3_1st <- model3_1st %>%
  mutate(
    coef_stars = paste0(round(coef, 3), add_stars(pval)),
    se_paren = paste0("(", round(se, 3), ")")
  ) %>%
  pivot_longer(cols = c(coef_stars, se_paren), names_to = "type", values_to = "estimate") %>%
  arrange(term, type) %>%
  dplyr::select(term, estimate)

model3_2nd <- model3_2nd %>%
  mutate(
    coef_stars = paste0(round(coef, 3), add_stars(pval)),
    se_paren = paste0("(", round(se, 3), ")")
  ) %>%
  pivot_longer(cols = c(coef_stars, se_paren), names_to = "type", values_to = "estimate") %>%
  arrange(term, type) %>%
  dplyr::select(term, estimate)

model_second_stage_m3 <- cbind(model3_1st, model3_2nd[,2])

modelsummary::datasummary_df(data = model_second_stage_m3, 
                             output = "../results/Table_A4_M2_2nd.docx", 
                             title = "Model 2", longtable = T)


## M1 Selection Model ##

sumamry_selection_m3 <- summary(selection_m3)

model4_1st_coef <- sumamry_selection_m3$tableP1[,1]
model4_2nd_coef <- sumamry_selection_m3$tableP2[,1]

model4_1st_se <- sumamry_selection_m3$tableP1[,2]
model4_2nd_se <- sumamry_selection_m3$tableP2[,2]

model4_1st_p <- sumamry_selection_m3$tableP1[,4]
model4_2nd_p <- sumamry_selection_m3$tableP2[,4]

# Model 2: First Stage 
model4_1st_coef <- bind_rows(lapply(model4_1st_coef, as.data.frame), .id = "model")
model4_1st_se <- bind_rows(lapply(model4_1st_se, as.data.frame), .id = "model")
model4_1st_p <- bind_rows(lapply(model4_1st_p, as.data.frame), .id = "model")

names(model4_1st_coef) <- c( "term", "coef")
names(model4_1st_se) <- c( "term", "se")
names(model4_1st_p) <- c( "term", "pval")

model4_1st <- model4_1st_coef %>%
  full_join(model4_1st_se, by = c("term")) %>%
  full_join(model4_1st_p, by = c("term"))

# Model 2: Second Stage
model4_2nd_coef <- bind_rows(lapply(model4_2nd_coef, as.data.frame), .id = "model")
model4_2nd_se <- bind_rows(lapply(model4_2nd_se, as.data.frame), .id = "model")
model4_2nd_p <- bind_rows(lapply(model4_2nd_p, as.data.frame), .id = "model")

names(model4_2nd_coef) <- c( "term", "coef")
names(model4_2nd_se) <- c( "term", "se")
names(model4_2nd_p) <- c( "term", "pval")

model4_2nd <- model4_2nd_coef %>%
  full_join(model4_2nd_se, by = c("term")) %>%
  full_join(model4_2nd_p, by = c("term"))

# Function to add stars based on p-values
add_stars <- function(pvalue) {
  ifelse(pvalue < 0.01, "***", 
         ifelse(pvalue < 0.05, "**", 
                ifelse(pvalue < 0.1, "*", "")))
}

model4_1st <- model4_1st %>%
  mutate(
    coef_stars = paste0(round(coef, 3), add_stars(pval)),
    se_paren = paste0("(", round(se, 3), ")")
  ) %>%
  pivot_longer(cols = c(coef_stars, se_paren), names_to = "type", values_to = "estimate") %>%
  arrange(term, type) %>%
  dplyr::select(term, estimate)

model4_2nd <- model4_2nd %>%
  mutate(
    coef_stars = paste0(round(coef, 3), add_stars(pval)),
    se_paren = paste0("(", round(se, 3), ")")
  ) %>%
  pivot_longer(cols = c(coef_stars, se_paren), names_to = "type", values_to = "estimate") %>%
  arrange(term, type) %>%
  dplyr::select(term, estimate)

model_second_stage_m4 <- cbind(model4_1st, model4_2nd[,2])

modelsummary::datasummary_df(data = model_second_stage_m4, 
                             output = "../results/Table_A4_M2_3rd.docx", 
                             title = "Model 2", longtable = T)

## Information Table ##

summary(selection_m1)
AIC(selection_m1)
BIC(selection_m1)
logLik(selection_m1)

summary(selection_m2)
AIC(selection_m2)
BIC(selection_m2)
logLik(selection_m2)

summary(selection_m3)
AIC(selection_m3)
BIC(selection_m3)
logLik(selection_m3)

#### Figure 6 #####

# extract summary information and coefficients
coefs_m1_outcome <- as.data.frame(summary(selection_m1)[["tableP2"]]) %>% 
  tibble::rownames_to_column("variable") %>% 
  dplyr::select(variable, est = "Estimate", std_err = "Std. Error", z_val = "z value", p_val = "Pr(>|z|)") %>% 
  as_tibble()

#Extract variance covariance matrix

matrix1 <- selection_m1$He[c(14:26),c(14:26)]

vcov_m1_outcome <- MASS::ginv(matrix1)

colnames(vcov_m1_outcome) <- rownames(matrix1)
rownames(vcov_m1_outcome) <- colnames(matrix1)

# Prepare simulations

sim_marginal1 <- function(coefs, vcov, temp.data1){
  bootsamples1 <- MASS::mvrnorm(1000, coefs$est, vcov)
  dim(bootsamples1)
  latent1 <- 1:nrow(temp.data1)
  for(i in 1:dim(bootsamples1)[1]){
    print(i)
    latent1 <- cbind(latent1,
                     (VGAM::probit(bootsamples1[i,pull(coefs[1,1])]*temp.data1[,"X1"]+
                                     bootsamples1[i,pull(coefs[2,1])]*temp.data1[,"X2"]+
                                     bootsamples1[i,pull(coefs[3,1])]*temp.data1[,"X3"]+
                                     bootsamples1[i,pull(coefs[4,1])]*temp.data1[,"X4"]+
                                     bootsamples1[i,pull(coefs[5,1])]*temp.data1[,"X5"]+
                                     bootsamples1[i,pull(coefs[6,1])]*temp.data1[,"X6"]+
                                     bootsamples1[i,pull(coefs[7,1])]*temp.data1[,"X7"]+
                                     bootsamples1[i,pull(coefs[8,1])]*temp.data1[,"X8"]+
                                     bootsamples1[i,pull(coefs[9,1])]*temp.data1[,"X9"]+
                                     bootsamples1[i,pull(coefs[10,1])]*temp.data1[,"X10"]+
                                     bootsamples1[i,pull(coefs[11,1])]*temp.data1[,"X11"]+
                                     bootsamples1[i,pull(coefs[12,1])]*temp.data1[,"X12"]+
                                     bootsamples1[i,pull(coefs[13,1])]*temp.data1[,"X13"], inverse=T)))
    
  }
  upper1 <- apply(latent1[,2:ncol(latent1)],1,quantile,.975)
  lower1 <- apply(latent1[,2:ncol(latent1)],1,quantile,.025)
  mean1 <- apply(latent1[,2:ncol(latent1)],1,mean)
  
  sims <- list(mean1,lower1,upper1,temp.data1)
  return(sims)
}

## Model 1 ##
# create data frame with coefficient values --> ResCap_lag
temp.data.ResCap_lag_m1 <- data.frame(`1` = 1,
                                  `2`= seq(min(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$ResCap_lag), max(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$ResCap_lag), by = 0.005),
                                  `3`= mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_gdppclog_lag, na.rm = T),
                                  `4` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_gdppc_growth_roll5_lag, na.rm = T),
                                  `5` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_poplog_lag, na.rm = T),
                                  `6` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$regionpol_EDI_lag, na.rm = T),
                                  `7` = 0,
                                  `8` = 0,
                                  `9` = 0,
                                  `10` = 0,
                                  `11` = 0,
                                  `12` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$year_1900, na.rm = T),
                                  `13` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$year_sq, na.rm = T))

##Model 2 ##

# extract summary information and coefficients
coefs_m2_outcome <- as.data.frame(summary(selection_m2)[["tableP2"]]) %>% 
  tibble::rownames_to_column("variable") %>% 
  dplyr::select(variable, est = "Estimate", std_err = "Std. Error", z_val = "z value", p_val = "Pr(>|z|)") %>% 
  as_tibble()

#Extract variance covariance matrix

matrix1 <- selection_m2$He[c(14:26),c(14:26)]

vcov_m2_outcome <- MASS::ginv(matrix1)

colnames(vcov_m2_outcome) <- rownames(matrix1)
rownames(vcov_m2_outcome) <- colnames(matrix1)

# create data frame with coefficient values --> MResCap_lag
temp.data.ResCap_lag_m2 <- data.frame(`1` = 1,
                                   `2`= seq(min(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$MResCap_lag), max(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$MResCap_lag), by = 0.005),
                                   `3`= mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_gdppclog_lag, na.rm = T),
                                   `4` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_gdppc_growth_roll5_lag, na.rm = T),
                                   `5` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_poplog_lag, na.rm = T),
                                   `6` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$regionpol_EDI_lag, na.rm = T),
                                   `7` = 0,
                                   `8` = 0,
                                   `9` = 0,
                                   `10` = 0,
                                   `11` = 0,
                                   `12` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$year_1900, na.rm = T),
                                   `13` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$year_sq, na.rm = T))


#### Model 3 #####

# extract summary information and coefficients
coefs_m3_outcome <- as.data.frame(summary(selection_m3)[["tableP2"]]) %>% 
  tibble::rownames_to_column("variable") %>% 
  dplyr::select(variable, est = "Estimate", std_err = "Std. Error", z_val = "z value", p_val = "Pr(>|z|)") %>% 
  as_tibble()

#Extract variance covariance matrix

matrix1 <- selection_m3$He[c(14:26),c(14:26)]

vcov_m3_outcome <- MASS::ginv(matrix1)

colnames(vcov_m3_outcome) <- rownames(matrix1)
rownames(vcov_m3_outcome) <- colnames(matrix1)


## Model 3 ##
# create data frame with coefficient values --> ResCap_lag
temp.data.ResCap_lag_m3 <- data.frame(`1` = 1,
                                      `2`= seq(min(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$AResCap_lag), max(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$AResCap_lag), by = 0.005),
                                      `3`= mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_gdppclog_lag, na.rm = T),
                                      `4` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_gdppc_growth_roll5_lag, na.rm = T),
                                      `5` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$e_poplog_lag, na.rm = T),
                                      `6` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$regionpol_EDI_lag, na.rm = T),
                                      `7` = 0,
                                      `8` = 0,
                                      `9` = 0,
                                      `10` = 0,
                                      `11` = 0,
                                      `12` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$year_1900, na.rm = T),
                                      `13` = mean(vdem_subset_v14_noNA_stage2[vdem_subset_v14_noNA_stage2$aut_ep==1,]$year_sq, na.rm = T))



# Simulating predicted probability #

sims_ResCap_lag_outcome_m1 <- sim_marginal1(coefs = coefs_m1_outcome, vcov = vcov_m1_outcome, temp.data1 = temp.data.ResCap_lag_m1)
sims_ResCap_lag_outcome_m2 <- sim_marginal1(coefs = coefs_m2_outcome, vcov = vcov_m2_outcome, temp.data1 = temp.data.ResCap_lag_m2)
sims_ResCap_lag_outcome_m3 <- sim_marginal1(coefs = coefs_m3_outcome, vcov = vcov_m3_outcome, temp.data1 = temp.data.ResCap_lag_m3)

data_sim_1 <- cbind(sims_ResCap_lag_outcome_m1[[4]][["X2"]], sims_ResCap_lag_outcome_m1[[1]])

# Plot effects for ResCap_lag
ResCap_lag_outcome_m1 <- ggplot() +
  geom_line(aes(x = sims_ResCap_lag_outcome_m1[[4]][["X2"]], y = sims_ResCap_lag_outcome_m1[[1]])) +
  geom_ribbon(aes(x =sims_ResCap_lag_outcome_m1[[4]][["X2"]], ymin = sims_ResCap_lag_outcome_m1[[2]], ymax = sims_ResCap_lag_outcome_m1[[3]]), alpha = 0.3) +
  ylab("P(breakdown resilience | autocratization)") + xlab("Resilience Capacity") + theme_bw() + 
  scale_y_continuous(limits = c(0.0,0.2), breaks = seq(0,0.5,0.1)) +
  scale_x_continuous(limits = c(0,1), breaks = seq(0,1,0.1)) +
  geom_rug(data = vdem_subset_v14_noNA_stage2, aes(x = ResCap_lag, y = 0), alpha = .25, position = "jitter", sides = "b")

df_outcome_m1 <- data.frame(x = sims_ResCap_lag_outcome_m1[[4]][["X2"]],
                            y = sims_ResCap_lag_outcome_m1[[1]])

# M2 
ResCap_lag_outcome_m2 <- ggplot() +
  geom_line(aes(x = sims_ResCap_lag_outcome_m2[[4]][["X2"]], y = sims_ResCap_lag_outcome_m2[[1]])) +
  geom_ribbon(aes(x =sims_ResCap_lag_outcome_m2[[4]][["X2"]], ymin = sims_ResCap_lag_outcome_m2[[2]], ymax = sims_ResCap_lag_outcome_m2[[3]]), alpha = 0.3) +
  ylab("P(breakdown resilience | autocratization)") + xlab("Multiplicative Resilience Capacity") + theme_bw() +
  scale_y_continuous(limits = c(0.0,0.2), breaks = seq(0,0.5,0.1)) +
  scale_x_continuous(limits = c(0,1), breaks = seq(0,1,0.1)) +
  geom_rug(data = vdem_subset_v14_noNA_stage2, aes(x = MResCap_lag, y = 0), alpha = .25, position = "jitter", sides = "b")



# Plot effects for ResCap_lag
ResCap_lag_outcome_m3 <- ggplot() +
  geom_line(aes(x = sims_ResCap_lag_outcome_m3[[4]][["X2"]], y = sims_ResCap_lag_outcome_m3[[1]])) +
  geom_ribbon(aes(x =sims_ResCap_lag_outcome_m3[[4]][["X2"]], ymin = sims_ResCap_lag_outcome_m3[[2]], ymax = sims_ResCap_lag_outcome_m3[[3]]), alpha = 0.3) +
  ylab("P(breakdown resilience | autocratization)") + xlab("Additive Resilience Capacity") + theme_bw() + 
  scale_y_continuous(limits = c(0.0,0.2), breaks = seq(0,0.5,0.1)) +
  scale_x_continuous(limits = c(0,1), breaks = seq(0,1,0.1)) +
  geom_rug(data = vdem_subset_v14_noNA_stage2, aes(x = AResCap_lag, y = 0), alpha = .25, position = "jitter", sides = "b")



### Resilience Capacity for Breakdown resilience ### check :)

figure <- egg::ggarrange(ResCap_lag_outcome_m1, ResCap_lag_outcome_m2, ResCap_lag_outcome_m3, nrow = 2, ncol = 2)

ggsave("../results/Figure_6.png", figure, width = 30, height = 30, units = "cm", dpi = 900, device = "png")


################################################################################
################################################################################

#### Figure 8 ####

#### Democratic Turnarounds ####

## lags for independent variables ##

vdem_subset_v14_noNA <- vdem_subset_v14 %>%
  mutate(e_gdppclog = log(e_gdppc+1), 
         e_poplog = log(e_pop +1), 
         year_1900 = year -1900) %>%
  group_by(country_id) %>%
  mutate(v2cademmob_lag = lag(v2cademmob), 
         e_gdppclog_lag = lag(e_gdppclog), 
         e_gdppc_growth_roll5_lag = lag(e_gdppc_growth_roll5), 
         regionpol_EDI_lag = lag(regionpol_EDI), 
         e_poplog_lag = lag(e_poplog), 
         v2x_polyarchy_lag = lag(v2x_polyarchy), 
         v2x_polyarchy_lag_squared = v2x_polyarchy_lag^2, 
         v2x_jucon_lag = lag(v2x_jucon), 
         v2xlg_legcon_lag = lag(v2xlg_legcon), 
         democracy_stock_lag = lag(democracy_stock), 
         ResCap_lag = lag(ResCap), 
         AResCap_lag = lag(AResCap), 
         MResCap_lag = lag(MResCap))

vdem_subset_v14_noNA <- vdem_subset_v14_noNA %>%
  mutate(region_factor = as.factor(e_regionpol_6C),
         region_factor = relevel(region_factor, ref=1), 
         spell = spell/100,
         spell_sq = spell^2/100,
         spell_cb = spell^3/100,
         year_sq = year_1900^2)


## drop NA ## 
vdem_subset_v14_noNA <- vdem_subset_v14_noNA %>%
  group_by(country_text_id) %>%
  mutate(v2x_regime_lag = lag(v2x_regime)) %>%
  ungroup() %>%
  filter(year >=2000) %>%
  mutate(aut_ep_stage2 = ifelse(aut_ep_outcome %in% c(1,5), 1, 0)) %>%
  drop_na(ResCap_lag, e_gdppclog, e_gdppc_growth_roll5_lag, regionpol_EDI_lag, 
          e_poplog,region_factor ,v2x_jucon_lag, v2xlg_legcon_lag, year_1900, spell, country_id, democracy_stock_lag) 


## List of all Autocratization Episodes ##

vdem_subset_v14_noNA <- vdem_subset_v14_noNA %>%
  group_by(country_text_id) %>%
  mutate(v2x_polyarchy_lag = lag(v2x_polyarchy)) %>%
  mutate(aut_ep_outcome_agg = case_when(aut_ep_outcome_agg == 1 ~ "Democratic breakdown", 
                                        aut_ep_outcome_agg == 2 ~ "No democratic breakdown", 
                                        aut_ep_outcome_agg == 3 ~ "Regressed autocracy", 
                                        aut_ep_outcome_agg == 4 ~ "Outcome censored")) %>%
  mutate(aut_ep_turnaround_id = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_id,
                                       ifelse(!is.na(uturn_ep_id), uturn_ep_id, NA)), 
         aut_ep_turnaround_start_year = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_start_year,
                                               ifelse(!is.na(uturn_ep_id), uturn_ep_start, NA)), 
         aut_ep_turnaround_end_year = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_end_year,
                                             ifelse(!is.na(uturn_ep_id), uturn_ep_end, NA)), 
         aut_ep_turnaround_outcome = ifelse(is.na(uturn_ep_id) & !is.na(aut_ep_id), aut_ep_outcome_agg,
                                            ifelse(!is.na(uturn_ep_id), uturn_ep_outcome, NA)))



figure_data_turnarounds <- vdem_subset_v14_noNA %>%
  group_by(country_text_id) %>%
  mutate(uturn_case = ifelse(any(aut_ep_turnaround_outcome == "J-shaped turnaround" | aut_ep_turnaround_outcome == "U-shaped turnaround" |
                                   aut_ep_turnaround_outcome == "L-shaped turnaround"), 1, 0)) %>%
  mutate(democracy_any = ifelse(any(democracy ==1), 1, 0)) %>%
  filter(uturn_case == 1) %>%
  filter(democracy_any == 1)

figure_data_turnarounds <- figure_data_turnarounds %>%
  mutate(aut_ep_turnaround_outcome = ifelse(is.na(aut_ep_turnaround_outcome), "no episode", aut_ep_turnaround_outcome))

ggplot(figure_data_turnarounds, aes(x= year, y = v2x_polyarchy)) +
  geom_line() +
  geom_line(data = subset(figure_data_turnarounds, aut_ep_turnaround_outcome == "J-shaped turnaround"), 
            aes(color = "J-shaped turnaround")) +
  geom_line(data = subset(figure_data_turnarounds, aut_ep_turnaround_outcome == "U-shaped turnaround"), 
            aes(color = "U-shaped turnaround")) +
  geom_line(data = subset(figure_data_turnarounds, aut_ep_turnaround_outcome == "L-shaped turnaround"), 
            aes(color = "L-shaped turnaround")) +
  geom_ribbon(aes(ymin = v2x_polyarchy_codelow, ymax = v2x_polyarchy_codehigh), alpha = 0.2) +
  geom_line(aes(y= ResCap), linetype = "dashed", color = "grey20") +
  geom_ribbon(aes(ymin = ResCap_codelow, ymax = ResCap_codehigh), fill = "grey20", alpha = 0.2) +
  facet_wrap(~country_name) +
  labs(color = "", 
       y = "Electoral Democracy Index/ \nDemocratic Resilience Capacity") + 
  theme_pubr() +
  theme(axis.text.x = element_text(size=10, angle=45)) +
  scale_color_manual(values = c("J-shaped turnaround" = "#00AFBB", "U-shaped turnaround" = "#E7B800", "L-shaped turnaround" = "#FC4E07")) +
  theme(legend.position="bottom") 


ggsave("../results/Figure_8.png", width = 25, height = 20, units = "cm", dpi = 900, device = "png")


#### Mean Resilience Capacity ####

mean_test_data <- vdem_subset_v14_noNA_stage2 %>%
  mutate(com_var = ifelse(aut_ep_turnaround_outcome == "J-shaped turnaround" & dem_ep ==1, "Democratic Turnaround",
                          ifelse(aut_ep_turnaround_outcome == "U-shaped turnaround" & dem_ep ==1, "Democratic Turnaround",
                                 ifelse(aut_ep_turnaround_outcome == "L-shaped turnaround" & dem_ep ==1, "Democratic Turnaround",
                                        ifelse(aut_ep_turnaround_outcome == "Democratic breakdown" & dem_breakdown_v2x_regime ==1, "Democratic Breakdown", NA)))))

mean_test_data <- mean_test_data %>%
  dplyr::select(country_text_id, country_name, year, com_var, everything())

table(mean_test_data$com_var)
mean_test_data <- mean_test_data %>%
  filter(!is.na(com_var)) %>%
  group_by(country_text_id, aut_ep_turnaround_id, com_var) %>%
  summarize(ResCap_lag = first(ResCap), 
            year = first(year))

## 43 countries as comparison

table(mean_test_data$com_var) ## Ten democratic breakdowns, 33 democratic turnarounds

ggplot(mean_test_data, aes(x = com_var, y = ResCap_lag, color = com_var)) +
  geom_boxplot(width=0.3) +
  geom_jitter() +
  labs(y = "Democratic Resilience Capacity",
       color = "", x = "",
       size = "",
       fill = "",
       shape = "",
       title = "",
       subtitle = "",
       caption = "") +
  theme_pubr() + 
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) +
  scale_color_manual(values=c('#8b3626', '#1F0561'))

ggsave("../results/Figure_7.png", width = 15, height = 10, units = "cm", dpi = 900, device = "png")





